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Abstract 

The  paper  describes  an  ongoing  work  with 
populating  the  world’s  largest  stress  intensity  factor  data 
base  with  92.4  million  new  solutions  and  separate  work 
consisting  of  large-scale  residual  strength  analysis  of  the 
C-130  Center-Wing-Box  (CWB)  considering  numerous 
different  multiple-crack  crack  configurations.  A 
computationally  efficient  and  reliable  procedure  is  used 
for  calculating  stress  intensity  factor  solutions  K(y)  to  be 
stored  in  the  data  base.  An  extended  technique  is  used  in 
predicting  the  residual  strength  of  the  C-130  CWB  for 
multiple  crack  configurations.  The  proposed  method 
requires  a method/solver  that  can  solve  the  very  complex 
nonlinear  contact  problems  between  rivets  and  the 
skin/stiffeners,  failure  of  rivets  with  a very  low 
computational  cost  per  crack  configuration.  The  splitting 
scheme  described  in  the  paper  is  the  basic  tool  used  to 
obtain  this  objective.  All  mathematical  equations  are 
solved  with  high  accuracy  with  respect  to  the  exact 
mathematical  solution  of  the  problem  and  with  control  of 
the  point-wise  error  (less  than  1%)  in  all  stress  intensity 
functions  K(y).  For  residual  strength  analysis  of  the  CWB, 
the  software  used  scales  very  well  on  computer  hardware 
like  SGI  Altix  (ASC/Eagle/Hawk)  and  IBM  P5 
(NA  VO/Babbage/Kraken) . Several  three-dimensional 
analyses  representative  of  the  size  and  complexity  of  the 
C-130  center  wing  box  have  been  completed.  An  example 
of  such  an  analysis  explicitly  modeled  the  wing  skins,  spar 
caps,  spar  webs,  and  stringers  which  resulted  in  90 
million  nodes  and  14  million  finite  elements.  Depending 
on  the  polynomial  order,  p,  used  in  the  solution,  the  total 
degrees  of  freedom  ranges  from  243-742  million  for 
polynomial  orders  p = 2 - 4;  respectively.  For  a more 
accurate  solution  polynomial  order,  5 is  needed  which 
results  in  a problem  with  1.2  billions  of  degrees-of- 
freedom. 


1.  Introduction 

State-of-the-art  fatigue  life  prediction  algorithms  can 
only  consider  rather  simple  structural  cracking  problems 
because  numerical  data  for  complex  multiple  crack 
configurations  are  simply  not  available.  Unfortunately 
modem  aircraft  are  complex  assemblies  with  diverse 
materials  and  joining  methods.  As  a result,  mission 
planners  and  combat  leaders  are  forced  to  use  much 
simplified  approaches  which  results  in  over-  or  un- 
conservative predictions  which  might  lead  to  costly 
decisions  and/or  jeopardize  aircraft  safety.  By  use  of 
novel  mathematical-numerical  methods  and  the 
Department  of  Defense  (DoD)  High  Performance 
Computing  Modernization  Program  (HPCMP) 
computational  resources,  the  situation  can  be  drastically 
improved  resulting  in  better  fleet  management  in  peace- 
and  war-time.  In  the  present  work  the  world  largest  stress 
intensity  factor  database  is  being  populated  with  high- 
accuracy  solutions.  Old  low  accuracy  solutions  are  also 
identified  and  replaced.  Today,  the  database  size  is  about 
two  orders  larger  than  that  available  in  2002  for  fatigue 
crack  growth  predictions.  This  is  a continuing  work.  The 
data  base  created  applies  to  both  civil  and  military 
aircraft.  Methods  and  data  developed  will  allow  the 
mission  planner  to  assess  the  level  of  risk  associated  with 
future  missions  with  respect  to  how  damage  accumulates 
per  flight.  The  present  work  also  covers  residual  strength 
analysis  of  the  CWB  of  the  C-130  aircraft  for  in-service 
observed  crack  configurations.  This  very  challenging  task 
is  briefly  summarized  below. 

The  present  paper  gives  details  with  regard  to 
creation  of  the  stress  intensity  factor  data  base.  In 
addition  we  describe  the  enabling  numerical- 
mathematical  technologies  for  residual  strength  analysis 
of  a full  C-130  CWB  where  all  geometrical  details  are 
modeled  as  three-dimensional  (3D)  objects  down  to 
rivet/bolt  size  and  nonlinear  behavior  captured  with 
control  of  error  in  the  numerical  solution. 
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2.  Worlds  Largest  K Database 

The  work  consists  of  deriving  92.4  million  stress 
intensity  factor  solutions  for  two  cracks  located  at  two 
holes  in  a plate  (Figure  1)  for  various  RJt,  a/t,  a/c,  D1/D2, 
L,  and  B,  for  through  thickness  cracks  and  not-through 
thickness  cracks  at  different  crack  locations. 

The  splitting  scheme  described  in  the  Users  Group 
Conference  2006  lecture1'1  was  adjusted  to  the  case  of  two 
cracks  with  crack  fronts  that  might  almost  be  touching. 
Figure  1 shows  schematically  the  three  different  solutions 
to  add  in  order  to  derive  the  solution  of  interest.  Each 
single  crack  problem  (totally  about  68,000)  is  analyzed 
separately.  From  the  single  crack  solutions  the,  92.4 
million  solutions  can  be  obtained  with  control  of  the  error 
as  described  in  Reference  1 and  references  therein.  The 
global  (crack  free)  problems  are  analyzed  for  specially 
designed  load  cases  where  traction  loading  and 
displacement  jumps  are  described  on  arbitrary  surfaces 
which  must  not  intersect  the  crack  surfaces.  Figure  2 
shows  one  such  surface  where  tractions  are  applied. 
Depending  on  the  crack  sizes,  different  loading  surfaces 
have  to  be  used  since  the  loading  surface  must  not 
intersect  a crack  surface.  Hence  each  one  of  the  68,000 
single  crack  problems  must  take  into  account  a hierarchy 
of  load  surfaces  which  depend  on  the  multiple  crack  cases 
of  interest. 

2.1.  Effectiveness  and  Robustness 

The  six-level  computational  scheme1'1  which  is  used 
to  derive  the  92. 4M  solutions  has  to  be  equipped  with 
automatic  functions  at  most  levels  if  manual  work  should 
not  be  prohibitively  large.  Examples  of  automatic 
functions  that  have  been  developed  are  mesh  generation 
(68K  meshes  for  the  /zp-version  of  finite  element  method 
[FEM]),  job  submission  (>68K  jobs),  job  monitoring 
functions,  scrubbing/resubmitting,  data  storage,  data 
transfer  etc. 

Submission  of  thousands  of  small  4-8  processor  jobs 
prevents  efficient  optimization  of  the  job  environment  to 
be  performed  by  the  system  load  scheduler.  Various  job 
clustering  techniques  have  therefore  been  developed  in 
order  to  support  system  optimization,  and  throughput  of 
the  68,000  jobs.  Significant  efforts  (i.e.,  direct 
simulation)  are  also  being  made  to  check  that  the  92.4M 
solutions  will  be  derived  with  a point  wise  error  less  than 
1%. 

3.  Residual  Strength  of  the'C-130  Hercules 
Aircraft  Center  Wing  Box 

Residual  strength  analysis  of  the  C-130  CWB 
requires  state-of-the-art  damage  and  fracture  analysis 


tools  combined  with  software  and  hardware  for  true  large- 
scale  analysis.  The  /zp-adaptive  code  STRIPE  which  has 
unique  capabilities  for  multi-scale  damage  and  fracture 
mechanics  analysis  is  used  for  this  purpose.  STRIPE  uses 
mixed  Open  Multi-Processing  (OpenMP)  and  Message 
Passing  Interface  (MPI)  and  executes  efficiently  on  the 
IBM  P4  system  KRAKEN/BABBAGE  at  the  Naval 
Oceanographic  Office  (NAVO)  and  at  the  SGI  Altix 
system  EAGLE  at  Aeronautical  Systems  Center  (ASC). 

The  main  questions  to  be  addressed  are  if  the 
C-130E/H  aircraft  now  being  retired  have  sufficient 
residual  strength  at  retirement  with  respect  to  fatigue 
cracking  in  the  CWB.  In  other  words,  was  the  aircraft  in 
service  too  long?  Conversely,  does  the  aircraft  have 
adequate  residual  strength  at  retirement  which  would 
allow  safe  usage  for  a significantly  longer  time;  i.e., 
aircraft  might  have  been  prematurely  retired. 

The  development  part  is  related  to  solution  of  very 
large  nonlinear  problems  of  the  size  of  billions  of 
degrees-of-freedom  (GDOF)  with  control  of  the 
numerical  error  in  the  solution.  Aircraft  structural 
mechanics  problems  of  this  complexity  have  hitherto 
never  been  solved.  Problems  of  such  complexity  remain  a 
target  in  industry  for  the  next  decade  in  order  to  do  more 
virtual  and  less  physical  testing. 

An  important  aspect  of  the  multi-scale  scheme  used 
is  efficient  handling  of  problems  with  several  100,000 
right  hand  sides1'1  for  problems  of  GDOFs  size. 
Significant  code  rewriting,  load  balancing  algorithms,  and 
code  optimization  have  so  far  resulted  in  an  efficiency 
which  for  problems  of  size  0.2  GDOFs  can  be  expressed 
as: 

Tlc  = T\  *(l+LC/J),  (1) 

In  Eq.  (1),  we  have  today /=4,000.  TLC(p ) is  the  wall 
time  needed  to  solve  the  problem  with  LC  right  hand  sides 
and  7)  the  wall  time  needed  to  solve  the  problem  with  one 
right  hand  side.  By  using  the  MPI-Input/Output  (IO) 
extensions  for  parallel  I/O  introduced  in  the  MPI-2 
standard  the  factor  / can  most  likely  be  significantly 
increased. 

Figure  3 summarizes  corrosion  damage  and  fatigue 
cracking  found  in  a retired  C-130  aircraft.  Cracks  were 
found  after  a complete  tear  down  process  where  all 
surfaces  were  inspected  for  cracks  and  corrosion  damage. 
Figure  3 shows  six  of  the  over  40  fatigue  critical  locations 
in  the  wing  box  and  also  regions  with  extensive 
corrosion121. 

Figure  4 shows  the  structural  domain  currently  being 
meshed  at  the  US  Air  Force  Academy  (USAFA).  A 
stringer  detail  is  shown  in  order  to  exemplify  the  mesh 
density  used.  The  mesh  generator  TrueGrid™,  developed 
by  XYZ  Scientific  is  used  to  create  the  large  finite 
element  mesh  which  is  estimated  to  consist  of  about  25 
million  finite  elements.  Derivation  of  accurate  solutions 
for  local  stresses,  including  control  of  the  error  in  the 
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solution  will  require  solution  of  problems  of  a size  1-2 
GDOF.  Structural  problems  of  this  size  are  about  two 
orders  larger  than  the  largest  structural  mechanics 
problems  being  solved  in  aircraft  industry  today. 

Efficient  solution  of  1-2  GDOF  problems  with 
several  100,000  right  hand  sides  constitutes  a true 
challenge.  Improved  domain  decomposition  division 
schemes,  load  balancing  algorithms,  memory/cache 
utilization  and  general  software  improvements  in  order  to 
handle  finite  element-meshes  consisting  of  about  25 
million  elements  has  been  developed. 

Several  generic  C-130  models  of  increasing  sizes 
were  created  and  successfully  analyzed  on  the 
ASC/EAGLE  system.  Figure  5 shows  a mesh  of  a 
generic  C-130-wing  consisting  of  7 million  hexahedral 
finite  elements.  A mesh  twice  that  large  was  solved  for 
polynomial  orders  p= 3 having  418  measured  degrees-of- 
ffeedom  (MDOF).  This  problem  might  be  the  largest 
aircraft  structural  mechanics  problems  ever  solved.  The 
work  continues  to  reach  objective  of  efficient  solution  of 
problems  of  GDOF-size.  Difficulties  encountered  are  in 
many  cases  believed  to  be  related  to  compiler  errors  and 
use  of  proper  environment  parameters  like  stack  size  limit 
etc. 

Table  1 summarizes  results  from  some  test 
computations  with  increasingly  larger  meshes  (only  64 
processors  were  used  in  the  tests).  Good  scaling  up  to 
512  processors  has  been  demonstrated. 

Table  1.  Benchmark  test  examples  executed  at 
ASC/EAGLE  (SGI  Altix)  using  64  processors  and  256 
GBYTE  of  memory 
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3.1.  Error  Estimation 

Despite  the  high  mesh  density  used,  the  error  in  local 
stress  has  somewhat  surprisingly  been  found  to  be  as 
large  as  15%  in  local  areas,  if  20  nodes  hexahedral  finite 
elements  are  used.  Figure  6 shows  a converged  nonlinear 
solution  for  p= 6 (105-noded  hexahedral  elements)  having 
a total  of  102  MDOF.  The  figure  shows  large  global 
buckles  interacting  with  local  buckles.  Near  local 
buckles,  red  areas  in  the  plot  in  Figure  6,  buckles  that 


might  contribute  to  de-bonding  between  stringer  and  skin 
and/or  rivet  failure  can  be  caused  by  the  15%  error  in  the 
p= 2 versus  p= 6 solutions..  The  /^-version  of  FEM  used 
in  the  present  project  makes  the  necessary  control  of  the 
error  in  the  solution  for  local  quantities  of  interest 
possible. 

3.2.  Accurate  and  Fast  Solution  of  3D  Contact 
Problems 

The  residual  strength  analysis  of  the  C-130  CWB 
requires  reliable  calculation  (i.e.,  error  control)  of  the 
maximum  load  carrying  capacity  of  the  wing  box  for  a 
given  multi-crack  scenario  with  interacting  cracks, 
buckling,  failing  rivets.  We  here  briefly  indicate  the 
strategy  used  to  solve  the  contact  problem  between 
rivets/bolts  and  the  aircraft  skin/stiffeners.  Due  to  the 
nonlinear  character  of  the  contact  problem  an  iterative 
scheme  has  been  designed  which  is  both  robust,  virtually 
exact,  converges  fast  and  simultaneously  takes  the 
mathematical  character  of  the  problems  appearing  during 
the  iteration  steps  into  account.  The  contact  problem  (see 
Figure  7 for  an  example  in  a two-dimensional  (2D) 
setting)  of  interest  is  split  into  three  problems,  similarly  as 
shown  in  Figure  1.  The  global  problems  a)  and  c)  in 
Figure  7 are  analyzed  without  considering  either  the 
contact  surfaces  or  the  crack  surfaces,  i.e.,  all  surfaces  are 
in  perfect  sliding-free  contact.  The  splitting  scheme  used 
then  has  the  advantage  that  the  global  analysis  a)  and  c) 
which  is  very  costly  for  problems  of  GDOF-size  do  not 
have  to  be  repeated  for  varying  crack  sizes  and  contact 
surfaces.  The  local  problems  b),  also  see  Figure  8,  are 
analyzed  for  fixed  (but  a priori  unknown)  locations  of  the 
contact  surfaces  for  a set  of  loading  functions  defined  by 
carefully  selected  traction  basis  functions  Q\  - Qt,  and 
scaling  factors  the  primary  unknowns. 

A necessary  condition  to  be  satisfied  at  the  two  points 
separating  contact  and  no-contact  is: 

K,=  0 (2) 

and  that  normal  stresses  are  compressive  at  the  part  of  the 
circular  boundary  being  in  contact.  Equation  (2)  can  be 
used  as  a very  precise  criterion  to  find  the  two  points 
separating  contact  and  no-contact.  A simple  Newton-type 
algorithm  can  then  be  used  to  find  the  contact  points 
satisfying  Eq.  (2)  for  each  bolt/rivet.  The  splitting 
scheme  is  simply  used  as  a very  fast  solver  for 
determination  of  stress  intensity  factors  at  the  crack-tips 
and  at  potentially  correct  contact  positions  during 
iterations.  Note  that  during  iterations  and  before 
convergence  is  achieved  Ki*  0 implying  that  there  is  a 
strong  stress  singularity  at  the  assumed  contact  positions. 
Robustness  of  the  approach  is  obtained  by  using  a hp- type 
of  mesh  moving  with  the  location  of  the  assumed  contact 
points.  Experience  shows  that  the  algorithm  converges  to 
about  three  digits  accuracy  in  less  than  five  iterations. 
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Figure  9 gives  an  overview  of  the  contact  solution  scheme  without  the  support  of  the  DoD  High  Performance 
in  a 3D  setting.  Computing  Modernization  Program. 

3.3.  Rivet  Failure  References 


Rivets  are  normally  not  loaded  in  tension  but  might 
see  tensile  loads  in  post-buckling  regime,  especially  in 
case  of  shear  buckling  but  also  in  presence  of  (larger) 
cracks.  The  model  used  to  simulate  bolt/rivet  failure  is 
based  on  a cohesive  zone  type  of  model  resting  on  an 
assumption  that  all  deformation  takes  place  in  a plane  in 
the  bolt.  Figure  10  exemplify  the  location  of  the  cohesive 
zone  in  case  of  a bolt.  The  adaptive  mesh  design  based 
on  the  ftp-version  of  FEM  which  follows  the  contact  lines 
are  also  indicated  in  Figure  10.  Normal  and  shear 
tractions  (er,  r)  are  assumed  to  depend  on  the  displacement 
jumps  over  the  cohesive  zone.  In  the  spirit  of  the  above 
mentioned  schemes  for  solution  of  crack  problems  and 
contact  problems  a modified  splitting  scheme  is  used  to 
simulate  rivet/bolt  failure  due  to  shear  and  tensile  loading. 
This  approach  maintains  the  robustness,  accuracy  and 
efficiency  of  the  overall  approach  and  makes  statistical 
analysis  economically  feasible. 

4.  Summary 

The  focus  of  the  present  paper  is  the  development  of 
advanced  computational  methods  for  reliable  fatigue  and 
residual  strength  analysis  of  bolted-riveted-bonded 
metallic  structures.  The  main  steps  discussed  are: 

• Creation  of  the  extended  data  base  where  this 
part  will  consist  of  92.4M  solutions 

• Fast  solver  for  solution  of  GDOF-problems  with 
several  100,000  right  hand  sides 

• Creation  of  FE-meshes  of  the  entire  C-130  main 
wing  box 

• Benchmarking  on  generic  C-130  models  with  up 
to  0.4  GDOFs 

• Development  of  a splitting  scheme  for 
simulation  of  rivet/bolt  failure  has  been  initiated 
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Figure  1.  Schematic  picture  of  splitting  scheme  in  a 20 
setting  and  in  case  stress  intensity  factor  calculation  for  the 
case  with  known  contact  surfaces  (see  Reference  1) 


Figure  2.  Definition  of  parameters  L,B,D1,D2  and  crack 
locations  1-4.  Two  cracks  at  locations  3 and  4 respectively 
are  shown  together  with  a proper  path  (blue  color)  for 
load/displacement  jump  application  for  this  crack 
combination. 
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• Fatigue  Critical  Locations:  CW 1,  CW 
5C,  CW  7B,  CW  9,  CW  1 4.  CW1 7 

•Corrosion  Indications:  SAIC  beam  cap 
corrosion  inspections 

•Crack  Indications:  SAIC  FWD  lower 
beam  cap  and  an  rainbow  fitting 
inspections  for  cracks 

Figure  3.  Fatigue  and  corrosion  sensitive  regions  found  after 
tear  down  and  inspection  of  a C-130  wing  box 


Figure  4.  Center  C-130  wing  box  which  currently  is  being 
meshed  in  large  detail  at  USAFA.  Stringer  shown  in  the 
figure  exemplifies  a characteristic  mesh  density. 


Figure  5.  Generic  '/i  mesh  of  C-130  wing  box  consisting  of  7 
million  hexahedral  elements 


Figure  6.  Interacting  local  and  global  buckles  in  generic 
fuselage.  The  converged  solution  shown  is  the  p=6  solution 
having  102  MDOF. 


Figure  7.  Schematic  picture  of  splitting  scheme  in  a 2D 
setting  and  in  case  stress  intensity  factor  calculation  for  the 
case  with  unknown  contact  surfaces 


Figure  8.  Loads  on  local  problems  are  shear  tractions  Qi  (<p) 
on  the  entire  circular  boundary,  normal  tractions  Q2  (<p)  on 
the  part  of  the  circular  boundary  not  in  contact,  shear 
tractions  Q3  ( tp ) and  normal  tractions  Q4  (<p)  on  the  crack 
face,  respectively.  Ail  loading  cases  are  treated  separately. 
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Figure  9.  Steps  in  solution  of  3D  contact  problems  using  a 
fracture  mechanics  approach  to  find  contact  areas  and  crack 
data  with  high  accuracy  in  a few  iterations 


Figure  10.  Cohesive  zone  model  used  for  robust,  accurate 
and  efficient  simulation  of  bolt/rivet  failure 
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